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ABSTRACT 

Context. Weak-lensing surveys need accurate theoretical predictions for interpretation of their results and cosmological- 
parameter estimation. 

Aims. We study the accuracy of various approximations to cosmic shear and weak galaxy-galaxy lensing and investigate 
effects of Born corrections and lens-lens coupling. 

Methods. We use ray-tracing through the Millennium Simulation, a large A'^-body simulation of cosmic structure for- 
mation, to calculate various cosmic-shear and galaxy-galaxy-lensing statistics. We compare the results from ray-tracing 
to semi- analytic predictions. 

Results, (i) We confirm that the first-order approximation (i.e. neglecting lensing effects beyond first order in density 
fluctuations) provides an excellent fit to cosmic-shear power spectra as long as the actual matter power spectrum is 
used as input. Common fitting formulae, however, strongly underestimate the cosmic-shear power spectra (by > 30% 
on scales £ > 10000). Halo models provide a better fit to cosmic shear-power spectra, but there are still noticeable 
deviations (~ 10%). (ii) Cosmic-shear B-modes, which are induced by Born corrections and lens-lens coupling, are 
at least three orders of magnitude smaller than cosmic-shear E-modes. Semi-analytic extensions to the first-order ap- 
proximation predict the right order of magnitude for the B-mode. Compared to the ray-tracing results, however, the 
semi-analytic predictions may differ by a factor two on small scales and also show a different scale dependence, (iii) The 
first-order approximation may under- or overestimate the galaxy-galaxy-lensing shear signal by several percent due to 
the neglect of magnification bias, which may lead to a correlation between the shear and the observed number density 
of lenses. 

Conclusions, (i) Current semi-analytic models need to be improved in order to match the degree of statistical accuracy 
expected for future weak-lensing surveys, (ii) Shear B-modes induced by corrections to the first-order approximation 
are not important for future cosmic-shear surveys, (iii) Magnification bias can be important for galaxy-galaxy-lensing 
surveys. 
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1. Introduction 

During the past few years, weak gravitational lensing has 
developed rapidl y from mere d etection to an important cos- 
mological tool ( Munshi et al.| [2008). Measurements of cos- 
mic shear help us to constrain tIie~properties of the cosmic 



matter distribution (e.g. |Semboloni et al. 


2006 Hoekstra 


et al. 


20061 ISimon 


et al.|2007||Benjamin et al.|2007||Massey 


et al. 


2007b 


'Fu e 


t al."2008'), the growth of structure (e.g. 
Massey et al. 2007c), and the nature of 


Bacon et al.j,2005 



Amendola et al.||2008j ). Weak galaxy-galaxy lensing can be 



used to study the properties of galactic dark-matter halos 
and the relation between luminous and dark matter (e.g. 
Mandelbaum et ar][2006, ,Simon et al.||2007] |Gavazzi et al.| 
2007t 



The accuracy that can be reached in weak-lensing sur- 
veys is determined by several factors. On the observational 



side, high accuracy requires large field sizes and deep obser- 
vations with a high number density of galaxies with mea- 
surable shapes. Moreover, it is crucial to obtain an accurate 
and unbiased measurement of galaxy ellipticities. Finally, 
for the interpretation of the resulting data and the in- 
ference of cosmological parameters, an accurate theoreti- 
cal model is needed. A thorough understanding of system- 
atic effects in weak lensing will become particularly im- 
portant with the advent of very large weak-lensing surveys 
such as CFHTL^ KID^ Pan-STARR^ and LSSTg) or 
the p lanned Dark Energy Survejj^ DUNE (Refregier et al. 
2006[ ), and SNAF0 For these surveys, the statistical uncer- 
tainties will be very small, so the accuracy will be limited 



^ http : / /www. cfht .hawaii . edu/Science/CFHLS 



http :// http: //www. astro-wise . org/projects/KIDS 



shilbertOastro . uni-bonn . de 



http : //pan-Starrs . if a. hawaii . edu 

http : / /www . Isst . org 
^ http : // www . darkenergy survey . org 
® http://snap.lbl.gov 
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by the remaining systematics in the data reduction and 
theoretical modehng. 

While significant improvement on image-e llipticity mea- 



surements are expected in the near future (Massey et al 



2007a), one still needs to investigate, how uncertain cur- 
rent theoretical predictions are, and how much improve- 
ment can be expected for these. Presently, the most ac- 
curate way to obtain predictions for weak-lensing surveys 
is to perform ray-tracing through large high-resolution A^- 
body simulations of cosmic structure formation (see, e.g. 



Wambsganss et al.|1998||Jain et al.|2000| [White fc Hu|2000 



Van Waerbeke et al. 2001! |Hamana fc M ellier'2001| [Vale fc| 
'White„2003, , White 2005^ . The drawback of this approach 
is that large TV- body simulations are computationally de- 
manding, so using them to explore the whole parameter 
space of cosmological models is currently unrealistic. On 
the other hand, ray-tracing simulations enable one to check 
the approximations and assumptions made in computation- 
ally less demanding (semi-) analytic models, and adjust and 
extend these models where necessary. 

Numerous ray-tracing methods have been developed to 
study the many aspects of gravitational lensing. Tree-based 



ray-tracing methods ( Aubert et aTp007 1 that adapt to the 
varying spatial resolution of A'-bociy simulations have been 
used to study the imp act of substructure o n strong lensing 



by dark matter halos (Peirani et al. 20081. Cluster strong 
lensing simulations, which require good mass modelling of 
galaxy clusters, usually ignore the matter distribution out- 



tracin 


g (e.g. Bartelmann & Weiss 


1994 


2007 


Rozo et al. 


2008p. 



1994 Meneghetti et al. 



Many simulations of weak lensing by clusters and large 



scale structure (e.g. [Wambsganss et al. 



2000 



Vale & White, 2003 



1998 Jain et al. 



Pace et al.| 2007p employ al- 
gorithms tha t are based on the mu ltiple-lens-plane ap- 



proximation (Blandford & Narayan 1986J to trace light 



Couchman et al. 


1999 Carbone et al. 20081 perform ray- 


tracing though t 
tial. In a simp] 


re three-dimensional gravitational poten- 
er approach (e.g. White & Vale 2004 


Heymans et al. ^ 


?006 Hilbert et al.„2007a|), the matter in 



paths onto a single lens plane, which is then used to calcu- 
late lensing observables. Recent simulations of CMB lensing 
use generalisations of the single- or multiple-plane approxi- 
mation that take the curvature of the sky into account (e.g . 
Das fc Bode|2008||Teyssier et al.|2008l|Fosalba et al.|2008| . 

In this work, we employ multiple-lens-plane ray-tracing 
through the Millennium Simulation (Springel et al.^ 2005, ) 
to study weak lensing. One of the largest A^-body simula- 
tions available today, the Millennium Simulation provides 
not only a much larger volume, but also a higher spatial 
and mass resolution than simulations used for earlier weak- 
lensing studies. In order to take full advantage of the large 
simulation volume and high resolution, the ray-tracing al- 
gorithm used here differs in several aspe cts from algorithms 
used in previous works (e.g. Jain et ar[2000 ). Here, we give 
a detailed description of our ray-tracing algorithm. 

Semi-analytic weak-lensing predictions are usually 
based on the first-order approximation, in which light de- 
flections are only considered to first order in the peculiar 
gravitational potential and hence, to first order in the mat- 
ter fluctuations. The ray-tracing approach allows us to look 
at effects neglected in the first-order approximation such 



Born corrections and lens-lens coupling. Here, we inves- 
tigate the cosmic-shear B-modes induced by these effects 



mates dCooray fc Hu 2002 


Hirata fc Seljak 2003 


Shapiro fc 


Coorayl 2006p, whose accuracy has not been confirmed by 



well fitti ng formu lae (Peac ock fc D odds 1996: Eisenstein &^ 
Hu 1999 ; [Smith et al. 20031 and halo models ( Seljak 2005[ 
Cooray fc Sheth 2002 i reproduce cosmic-shear power spec- 
tra. Finally, we investigate the accuracy of the first-order 
approximation for weak galaxy-galaxy lensing. 

The paper is organised as follows. In the next section, 
we introduce the theoretical background and notation used 
in our lensing analysis. In Sec. [3] we discuss our ray-tracing 
algorithm. The results from our ray-tracing analysis are 
presented in Sec.|4] We conclude our paper with a summary 
in Sec. |5] 

2. Theory 

2.1. Gravitational light deflection 

In this section, we introduce the formulae relating the 'ap- 
parent' positions of distant light sources to their 'true' posi- 
tions. In order to label spacetime points in a model universe 
with a weakly perturbed Friedmann-Lemaitre-Robertson- 
Walker (FLRW) metric, we choose a coordinate system 
(t, /3, w) based on physical time i, two angular coordinates 
(3 = (/3i, /32), and the line-of-sight comoving distance w rel- 
ative to the observer. The spacet ime metric of the model 
universe is then given by (see, e.g., Bartelmann fc Schneider 
20011): 



d.^ = fi+^uw-ri-^'j>\ 



c^ / \ c^ 

x\dw^ +fl{w) [dPl +cos2(/3i)d/32] \. (1) 



Here, c denotes the speed of light, a = a{t) denotes 
the scale factor, and <i> = ^{t,f3,w) denotes the peculiar 
(Newtonian) gravitational potential. The comoving angu- 
lar diameter distance is defined as: 



for K >0, 

for X = 0, and (2) 
for K <0, 



l/Viirsin 
fK{w) = { w 

,1/x/^sinh 



Kw 



where K denotes the curvature of space. The particular 
choice for the angular coordinates /3 = {^1,(^2) is con- 
venient for the application of the 'flat-sky' approxima- 
tion, where the metric near /3 = is approximated using 
cos2(/3i) « 1. 

Consider the path, parametrised by comoving distance 
ui, of a photon eventually reaching the observer from angu- 
lar direction 0. The angular position j3{6^w) of the photon 
at comoving distance w is then given by (see, e.g., Jain fc 
Seljak 1997, for a sketch of a derivation): 

P[0, w) = ^ / aw 



fK{w)jK{w') 

xVp<i>{t{w'),f3{e,w'),w') (3) 

with V/3 = ((9/9/3i, 9/9/32), and t{w') denoting the cosmic 
time of events at line-of-sight comoving distance w' from 
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the observer. By differentiation this equation w.r.t. 6, we 
obtain the distortion matrix A, i.e. the Jacobian of the lens 
mapping 6 t-^ (3 — /3(0, w): 



= Sii - 



fK{w - w') 



fK{w)fK{w') 



(4) 



Akj{e,w'). 



Due to the matrix products in Eq. Q, the distortion 
matrix A is generally not symmetric. However, it can be 
decomposed into a rotation matrix (related to a usually 
unobser vable rotation in the so urce plane) and a symmetric 
matrix ( Schneider et al.||l992 l: 



cos uj sm uj 
— sin u! cos UJ 



K - 71 

-72 1 



-72 
K + 7i 



(5) 



The decomposition defines the rotation angle uj = uj{0, w), 
the convergence k — k{6, w), and the two components 71 = 
ji{6,w) and 72 — 72(^1 w) of the shear, which may be 
combined into the complex shear 7 = 71 + 172. 

The shear field j{0, w) can be decomposed into a 
rotation-free part 7e(0, and a divergence-free part 
7b(0,u;). For infinite fields, the decomposition into these 
E/B-modes is most easily written down in Fourier space: 

7e(^, = [(^? - ^2)71 w) + 2^i^272(^, w)] , (6a) 

%{e,w) = ^ [{el-il)u^,w) - 2e,i2^,{£,w)] . (6b) 

Here, hats denote Fourier transforms, £ — {£1,12) denotes 
the Fourier wave vector, and £ = ii + i£2- Care must be 
taken when decomposing the shear in fields of fi nite size, 
where the field boundaries can cause artifacts (Seitz & 



Schneider 19961. These artifacts can be avoided by usmg 



aperture m asses to quantify the shear E- and B-mode co n- 
tributions (Critt enden et al.||2002[ [Schn eider et al.||2002[ ). 

Equations (|3| and Q are implicit relations for the light 
path and the Jacobian. The solution of Eq. ^ to first or- 
der in the potential is obtained by integrating along undis- 
turbed light paths: 

Jo fK{w)fK{w') 

xVg<i>{t{w'),e,w'). (7) 
The distortion to first order reads: 



2 / 

Aij{6,w) = Sij dw 



J fK{w - w') 
fK{w)fK{w') 

mW, • 

The first-order approximation to the distortion contains the 
Born approximation, which ignores deviations of the ac- 
tual light path from the undisturbed path on the r.h.s. of 



Eq. Moreover, lens-lens coupling is neglected, i.e. the 
appearance of the distortion on the r.h.s. of Eq. The 
neglected lens-lens coupling and corrections to the Born ap- 
proximation account for the effect that light from a distant 
source 'sees' a distorted image of the lower-redshift matter 
distribution due to higher-redshift matter inhomogeneities 
along the line-of-sight. Thus, the first-order approximation 
works well in regions where larger matter inhomogeneities 
are absent or confined to a small redshift range, but fails 
in regions where noticeable distortions arise from matter 
inhomogeneities at multiple redshifts. 

Born corrections and lens-lens coupling effects may cre- 
ate shear B-modes. The perturbative calculation of the 
she ar B-modes by iteratively solving Eg. Q is possi- 
ble ( [Cooray fc Hu|2002[ plirata fc Seljak|2003| ), but tedious, 
and the accuracy of this approach is not known. However, 
multiple deflections and lens-lens coupling effects are fully 
included in the multiple-lens-plane approximation as de- 
scribed below. We will thus use this approximation to in- 
vestigate these effects and assess the quality of perturbative 
calculations of these effects. 



2.2. The multiple-lens-plane approximation 

In the multiple-lens-plane approximation (see, e.g., 
Blandford fc Narayaril [19861 [Schneider et allpJ92| [Seitz 



et al. 1994 Jain et al. 2000), a series of lens planes per 
pendicular to the central line-of-sight is introduced into the 
observer's backward light cone. The continuous deflection 
that a light ray experiences while propagating through the 
matter inhomogeneities in the light cone is then approxi- 
mated by finite deflections at the lens planes. The deflec- 
tions are calculated from a projected matter distribution 
on the lens planes. This corresponds to solving the inte- 
gral equations (jsj) and Q by discretisation (and using the 
impulse approximation) . 

The deflection oS^\l3'^^^) of a light ray intersecting the 
fc"^ lens plane (here, we count from the observer to the 
source) at angular position /J^'^) can be expressed as the 
gradient of a lensing potential -(/;'■*'•': 



c(fc)(/3(fc)) ^ V^(fe) ^/.('^)(/3('=)) 



(9) 



The differential deflection is then given by higher deriva- 
tives of the lensing potential. The second derivatives can 
be combined into the shear matrix 



) _ a^^(^)(/3(^)) _ 5af^(/3W) 



dp. 



(fc) 



(10) 



The lensing potential tp'^^'^ is a solution of the Poisson equa- 
tion: 

V^(,)V'"H/3'") = 2aW(/3^'^)- (11) 

The dimensionsless surface mass density cr*^'^) is given by a 
projection of the matter distribution in a slice around lens 
plane: 



a 



2c2 



dw'6^{f3'^''\w'). (12) 



(fc) 



Here, Hq denotes the Hubble constant, flm the mean mat- 
ter density in terms of the critical density, = /k(w*-'^') 
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and a^*^^ — a(u'''^''), with w*^*^) denoting the hne-of-sight 
comoving distance of the plane. Furthermore, Sjn{f3^'^\ if') 
denotes the three-dimensional density contrast at comov- 
ing position (^l3^''\'w'^ relative to the mean matter den- 



sity. The slice boundaries 



,(fc) 



< w 



(k) 



< w 



(fc) 



and 



„(fc) 



and 



,(fe) 



have to satisfy 
They are usually 



chose n to correspond to the mean redshifts (e.g. Jain et al 



2000 1 or comoving distances (e.g. Wambsganss et al.||2004p 



of successive planes These conditions ensure that every re- 
gion of the light cone contributes exactly to one lens plane, 
which is the closest plane in redshift or comoving distance. 

Given the deflection angles on the lens planes, one can 
trace back a light ray reaching the observer from angular 
position (3''^^ ~ on the first lens plane to the other planes: 



fe-l r{i,k) 

/3W(0) = 0-^:^a«^ 



■Ik 



(M')), A; -1,2,., 



(13) 



Here, /^'"^ ^ fiAuM) 



Equation ( |13[ ) is not practical for tracing rays through 
many lens planes. An alternativ e e xpression is obtain ed as 
follows (see, e.g., Hartlap 2005 or Seitz et al. 1994 for a 



different derivation): The angular position P'^'^^ of a light 
ray on the lens plane k is related to its positions (3^''^'^^ 
and f)''-''^^^ on the two previous lens planes by (see Fig.[l]): 



where e 



+ ft 



2,k) 



Ak- 
Jk 



(14) 



Ak- 

J K 



Ak- 

J K 



Hence, 



/3 



(fe) 



Ak 



K 



-1) 



Ak- 
J K 



(k-1) 



/3 



-2,k) 



Ak) 
J K 



Ak- 

J K 



2,/c-l) 



/3 



(fc-2) 



IK 



Ak-2.k) 
■IK 



Ak) 
J K 

Ak- 

J K 

Ak) 
J K 



Ak- 

J K 



2, fe-l) 



/3 



(fe-i) 



(15) 



l,k) 



a 



For a light ray reaching the observer from angular position 
on the first lens plane, one can compute its an gula r po- 
sition on the other lens planes by iterating Eq. (15 1 with 
initial values /3(°) = /3(i) = 0. 

Differentiating Eq. ^TE\ with respect to 0, we obtain a 
recurrence relation for the distortion matrix: 



Ak) 



Ak 



K 



Ak 
■I i 



-2,fe) 



K 



Ak-2) 



Ak) 
J K 



Ak- 
J K 



2,k- 



Ak-1) 
J K 



Ak) 
J K 
Ak- 
J K 
Ak) 
J K 



Ak-2,k) 

[k /i,ik-l) 

1) 



Ak 

J K 



-2.k~ 



(16) 



l,k) 



.,(fc-l) .(fc-1) 
'^ik "fej 



^ The exact choice for the projection boundaries becomes 
unimportant for sufficiently small spacings between the lens 
planes. 



With the knowledge of the involved distances and shear 
matrices, this equation allows us to iteratively compute the 
distortion matrix of a light ray from the observer to any 
lens plane. This equation requires in practice much fewer 
arithmetic operat ions and memory than the commonly used 
relations (e.g. by Jain et aL||2000 ) based on Eq. (13 1. 

For comparison and testing, we will also use the 
multiple-lens-plane algorithm to calculate the distortion in 
the first-order approximation by: 



fe — 1 r{n.k) 

E JK ||(") 
Ak) ^ij 
n=l J K 



{6). 



(17) 



3. The ray-tracing algorithm 

The methods we use for ray-tracing through A'^-body simu- 
lations t o study len s ing ar e g enerally similar to tho se used 



by, e.g., Jain et al. (20001 or Vale & White (20031. First 



the matter distribution on the past light cone of a fiducial 
observer is constructed from the simulation data. Then, the 
past light cone is partitioned into a series of redshift slices. 
The content of each slice is projected onto a lens plane. 
Finally, the multiple-lens-plane approximation is used to 
trace back light rays from the observer through the series 
of lens planes to the sources. 

The purpose of our ray-tracing algorithm is to simulate 
strong and weak lensing in a way that takes full advan- 
tage of the unprecedented statistical power offered by the 
large volume and high spatial and mass resolution of the 
Millennium Simulation {^Therefore, our ray-tracing method 
differs in many details from previous works. Most notably, 
we use a multiple-mesh method and adaptive smoothing 
to calculate light deflections and distortions from the pro- 
jected matter distribution on the lens planes. This allows 
us to simulate lensing on the full range of scales covered by 
the Millennium Simulation, ranging from strong lensing on 
scales > 1 arcsec to cosmic shear on scales < 1 deg. A brief 
outline of our algorithms for the construction of the past 
light cones and the lens pl anes has been given in an earlier 
work ( Hilbert et al.|2007b l. Here, we extend the discussion 
and provide a more detailed description. 

3.1. The Millennium Simulation 



The Millennium Simulation (Springel et al. 2005) is a 



large A^-body simulation of cosmic structure formation in 
a flat ACDM universe. The following cosmological param- 
eters were assumed for the simulation: a matter density of 
r^ni = 0.25 in units of the critical density, a cosmological 
constant with J^a = 0.75, a Hubble constant h = 0.73 in 
units of 100kms~^Mpc~^, a spectral index n = 1 and a 
normalisation parameter erg = 0.9 for the primordial lin- 
ear density power spectr um. These chosen p arameters are 
consistant with the 2d F dCoUess et al.|l2003| ) and WMAP 

The Simula- 



Ist-year data analysis ( Spergel et al 

tion followed the evolution of the matter distribution in 
a cubic region oi L — 500/i^^Mpc comoving side length 
from redshift z — 1 27 to t he present using a TreePM ver- 
sion of GADGET-2 (ISpringel 2005) with 2160^ particles of 



This work concentrates on weak lensing, but the algorithm 



is also used for str ong-lensing studies ( Hilbert et al. 2007b 
Faure et al.|2008[ ) 



2008 
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1 ... k-2 k-l k 




rW r(k-2) (k-l) (k-l) (k-l) (k) 

Jk ••• Jk Jk,l Jk Jk,\] Jk 



Fig. 1. Schematic view of the observer's backward Hght cone in the multiple-lens-plane approximation. A light ray (red 
line) experiences a deflection only when passing through a lens plane (solid blue lines). The deflection angle a^''^^^ of a 

(k—l) 

ray passing through the lens plane at distance from the observer is obtained from the matter distribution between 

/if u^^ ^^'^ /if L^^ projected onto the plane. Using the deflection angle a^*-'"^) of the hght ray at the previous lens plane 
and the ray's angular positions /S^*^"^^ and /3('^~^) on the two previous planes, the angular position /S^*^) on the current 
plane can be computed. 



mass nip = 8.6 x 10 



h~ Mq and a force softening length 
of 5h~^ kpc comoving. 

Snapshots of the simulation were stored on disk at 64 
output times. These snapshots contain, among other data, 
the positions, SPH smoothing lengths and friend-of- friend 
(FoF) group data of the particles. The storage order for the 
particle data is based on a spatial oct-tree decomposition of 
the simulation cube (see Springel|2005 Springel et al.|2005" 



for details), which facilitates access to the particle data for 
small subvolumes of the simulation. 

Complex physical processes of baryonic matter such as 
the formation and evolution of stars in galaxies has not 
been incorporated directly into the Millennium Simulation. 
However, several galaxy- formation models have been used 
to predict the propertie s of galaxies in the simulation 
^Springel et al |[2b05 Crot on et al.| 2006 1 [Bower et al.|2006" 
Ue Lucia fc fj iaizot| |2007).THe" ray-tracing methods pre- 



sented in this paper will allow us (in future work) not only 
to study cosmic shear in great detail, but also to make pre- 
dictions for galaxy-galaxy lensing (and related higher-order 
statistics) for the various galaxy-formation models. 

3.2. The construction of the matter in the backward light 
cones 

Even with a comoving size of L = 500/i~^ Mpc, the simu- 
lation box is too small to trace back light rays within one 
box. We therefore exploit the periodic boundary conditions 
of the simulation by arranging replicas of the simulation 
box in a simple cubic lattice with a lattice constant equal 
to the box size L to fill space. We refrain from randomly 
shifting or rotating the content of the lattice cells, because 
the simulation box is far too large to be projected onto a 
single lens plane. In addition, this allows us to keep the 
matter distribution continuous across the cell boundaries. 

In this periodic matter distribution, light rays would en- 
counter the same structures many times at different epochs 
before reaching relevant source redshifts if one chose the 
line of sight (LOS) parallel to the box edges. Hence, the 



LOS must be chosen at a skewed angle relative to the box 
axes. On the other hand, the application of Fourier methods 
for the calculation of the light deflection at the lens planes 
requires a matter density that is periodic perpendicular to 
the LOS. Choosing a LOS parallel to n = (ni, n2, ns) with 
suitable coprime rij € Z, one can obtain a large enough 
repetition length of \n\L along the LOS (see Appendix [A| . 
At the same time, the matter distribution is periodic per- 
pendicular to the LOS with an area of periodicity given 
by |n|L^. Our choice for n = (1,3,10) yields a LOS pe- 
riodicity of 5.24/i~^ Gpc (corresponds to 2; = 3.87) and a 
rectangular unit cell of 1.58/i~"'^ Gpc x l.&Qh~^ Gpc for the 
lens planes. Moreover, any directions with shorter period- 
icity are at least 1.81 deg away from n, and a light cone 
with a 1.7 deg x 1.7 deg field of view does not intersect with 
itself up to redshift z — 3.87 when folded back into the 
simulation cube|^The resulting orientation of the LOS and 
the lens planes w.r.t. the simulation box are illustrated in 
Fig.|2] 

We partition the observer's backward light cone into 
redshift slices such that each slice contains the part of 
the light cone that is closer in redshift to one of the 
snapshots than any other snapshot (with exceptions near 
the boundaries discussed below). The boundary between 
two redshift slices with snapshot redshifts z^*^^ and z^'^"'"^) 



is thus a plane at comoving distance w 



{k) 
u 



„(fc+i) 



w [(zC') -I- zC^+i)) /2]. In addition, w"^' = 0. The particle 
data of the snapshot closest in redshift is then used to ap- 
proximate the matter dist ribution in each of the se slices. 
Fast box-intersection tests ( Gottschalk et al.|1996 l and the 
spatial-oct-tree storage order of the simulation are utilised 
to minimise reading of the particle data (which reduces run 
time by factors 5-10). 



^ We often use a larger field of view - in particular, if only 
lower source redshifts are considered. Even for high source red- 
shifts, where the resulting light cone may cover the same simu- 
lation region more than once, a lar ge field of view can be used 
with due care ( [Hilbert et al.|2007a[ ) 
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(a)z=Zi. : (b)z=Zi+i: 





Fig. 2. Schematic view of the orientation of the hne-of- 
sight (red hne) and the lens planes (blue area) relative to 
the simulation box (indicated by black lines). 



In the construction of the matter distribution in the 
light cone, special care is taken for the particles near the 
boundary of two slices. In the simulation, particle concen- 
trations representing dark matter halos of galaxies or clus- 
ters were identified by a friend-of-friend (FoF) group find- 
ing algorithm. Some of these halos are located on the slice 
boundaries with particles on either side |^ In order to avoid 
that such a halo is only partially included into a slice (and 
hence would be only partially projected onto a lens plane) , 
a halo is either included as a whole if its central particle 
is inside the slice as defined by boundary planes, or com- 
pletely excluded otherwise. 

If the matter structure in the simulation were static, 
this procedure would suffice to prevent parts of the same 
halo from being projected onto adjacent lens planes, which 
would create artificial close pairs of halos on the sky. Halos, 
however, may have moved across a slice boundary between 
two snapshots. We therefore amend the above inclusion cri- 
terion for halos near the boundary: If a halo is included in 
(excluded from) the slice of the later snapshot, its progen- 
itors in the earlier snapshot are excluded from (included 
in) the 'earlier' slice even if their centres lie on the 'early' 
('late') side of the slice boundary. These inclusion criteria 
for halos are illustrated in Fig. [3j 

3.3. The lens planes 

The matter content of each redshift slice of the backward 
light cone is projected along the LOS direction onto a lens 
plane. Each lens plane is placed at the comoving distance 
of the corresponding snapshot's redshift. The lens planes 
serve also as source planes for the ray-tracing. The resulting 
number of lens planes as a function of the source redshift 
is shown in Fig. |4j 

The light deflection angles and distortions resulting 
from the projected matter density on the len s planes are 
computed by p article-mesh (PM) methods (Hockney & 



Eastwood 
once the d 



1981). Mesh methods have the advantage that 
eflection and distortion are computed on a mesh 



Approx. 0.5% (5%) of halos with virial masses 
M200 > lO^^/i"^M0 (> IO^^Ii'^Mq) are affected by this pro- 
cedure. Though not essential for cosmic shear simulations (test 
show a relative difference ~ 0.1% for the shear power spectra), 
the proper treatment of halos near slice boundaries is important 
for group-galaxy lensing and strong lensing simulations. 




matter 
evolution 




■ in slice k 
I I in slice k+l 



Fig. 3. Schematic view of the adaptive slice boundaries to 
avoid the truncation or double inclusion of halos that are 
located near a slice boundary. Halos near the boundary of 
slice k and fc -t- 1 are either included as a whole in slice k 
or completely excluded depending on the positions of their 
centres (a). Halos that are included (excluded) in slice k, 
are excluded (included) from slice fc -I- 1 even if they have 
crossed the slice boundary between redshift k and k+l (b). 




Fig. 4. The number Nj^ of lens planes used for the ray- 
tracing as a function of the source redshift zs- 



(e.g. by Fast Fourier methods), the computation of the de- 
flections and distortions for many light rays intersecting 
the plane is very fast (compared to, e.g., direct-summation 
or tree methods). One disadvantage is that the used mesh 
spacing limits the spatial resolution of the projected mat- 
ter distribution. However, any A^-body simulation provid- 
ing the matter distribution for the ray-tracing has a limited 
resolution as well. In dense regions, the spatial resolution of 
the Millennium Simulation is effectively determined by the 
force softening, which is 5h~^ kpc comoving. Thus, a mesh 
spacing of 2.5h~^ kpc comoving is required to avoid resolu- 
tion degradation for the projected matter density. However, 
a single mesh covering the full periodic area of the lens 
plane (i.e. 1.58h^^ Gpc x 1.66/i^^ Gpc comoving) with such 
a small mesh spacing would be too demanding, in particu- 
lar regarding the memory required both for its computation 
and storage. We therefore use a hierarchy of meshes instead. 

The lensing potential '0 is split into a long-range part 
'0iong and a short-range part V'short- The split is defined in 
Fourier space by: 



V'short(^) = m [1 - exp 



and 



(18) 
(19) 



The splitting angle P^pUt = rspia/ fxiw), with comoving 
splitting length rgput and comoving angular diameter dis- 
tance of the lens plane fK{w), quantifies the spatial scale of 
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the spUt. Different meshes are then used to calculate V'long 
and -f/jshort- 

First, the particles in each slice are projected onto a 
coarse mesh of 16384 x 16384 points covering the whole 
periodic are a of the lens plane using clou ds-in-cells (CIC) 
assignment (Hockney & Eastwood 1981). The long-range 
potential V'long is then calculated on this mesh by means 
of fast Fou r ier transform (FFT) t echniques ( Cooley & 



Tukey||1965l |Frigo fc Johnson||2005[ ). The splitting length 
Tspiit = 175/1 kpc is chosen slightly larger than the coarse 
mesh spacing {96h^^ kpc and 101h~^ kpc comoving, respec- 
tively), so the coarse mesh samples "(Along with sufficient 
accuracy. For each lens plane, the long-range potential is 
calculated once, and the result is stored on disk for later 
use during the ray-tracing. 

The short-range potential i/'short is calculated 'on the 
fly', i.e. during the actual ray-tracing. The area where the 
light rays intersect the plane is determined and, if larger 
than 40/i~^ Mpc comoving, subdivided into several patches 
up to that size. Each patch is covered by a fine mesh 
with a mesh spacing of 2.5/i~^kpc comoving and up to 
16384 X 16384 mesh points. The fine meshes are chosen 
slightly larger than the patches in order to take into account 
all matter within the effective range of '(/'short, for which we 
assume 875/i~^kpc {— Srgpnt). The limited range of V'short 
ensures that the matter distribution outside the mesh af- 
fects only mesh points close to its boundary (i.e. within the 
effective range), but not the interior mesh points used for 
the subsequent analysis. Periodic boundary conditions can 
therefore be used for the FFT on the patches without 'zero 
padding'. 

In order to reduce the shot noise from the individual 
particles, either a fixed or an adaptive smoothing scheme is 
used for the matter distribution on the fine meshes. In case 
of the fixed smoothing, the particles in the slice are pro- 
jected onto the fine mesh using CIC. The resulting matter 
density on the fine mesh is then smoothed in Fourier space 
with a Gaussian low-pass filter 



Ks{£;Ps) = exp 



(20) 



whose filter scale /3s = Is/fxiw) is determined by the lens 
plane's comoving distance w and a fixed comoving filter 
scale Is- This is done during the calculation of the short- 
range potential "(Ashort with FFT methods. 

In case of the adaptive smoothing, the mass associated 
with each simulation particle contributes 



Sp(a;) 



3mn 



0, 



\X Xp\ < Tp, f^2l^ 
\x Xp I ^ Tp , 



to the surface mass density on the fine mesh. Here, x de- 
notes comoving position on the lens plane, Xp is the pro- 
jected comoving particle position, and Tp denotes the co- 
moving distance to the 64th nearest neighbour particle 
in three dimensions (i.e. before projection). The adaptive 
smoothing is essentially equivalent to the assumption that, 
in three-dimensional space, each simulation particle repre- 
sents a spherical cloud with a Gaussian density profile and 
an rms radius that is half the distance to its 64th nearest 
neighbour. From the resulting surface mass density on the 



fine mesh, the short-range potential ■(/'short is then calcu- 
lated by FFT methods. 

The long- and short-range contributions to the deflec- 
tion angles and shear matrices are calculated on the coarse 
and fine mesh by finite differencing of the potentials The 
values between mesh points are obtained by bilinear inter- 
polation. The resulting defiection angles and shear matrices 
at the ray positions are then used to advance the rays and 
their associated distortion matrices from one plane to the 
next. 

The ray-tracing algorithm reproduces the defiection an- 
gles and distortions caused by a single point mass very ac- 
curately, as is shown in Fig. [5] and |6] For the defiection 
angle, the relative deviation from the known analytical re- 
sult is at most one percent, with a much smaller rms error. 
Apart from scales below 10/i^^kpc comoving, where dis- 
creteness of the fine mesh becomes important, the largest 
relative errors occur on the scale of the split between short- 
and long-range potential. If desired, an even smaller error 
in this region could be obtained by increasing the splitting 
scale. 

In this work, we do not consider the effects of the stel- 
lar mass in galaxies. Note, however, that the ray-tracing 
algorithm can be extended to include the gravitational ef- 
fects of the stars in galaxies as described in Hilbert et al. 
(2008): The positions and stellar masses of the galaxies are 
inferred from semi-analytic galaxy-formation models imple- 
mented within the evolv ing dark-matter distribution of the 
Millennium Simulation ( Springel et al.||2005| Croton et al. 



2006) |De Lucia fc Blaizot||2007| ). I'he Tight deflection in 



duced by the stellar matter is then calculated using ana- 
lytic expressions for the projected stellar mass proflles on 
the lens planes. The error made by adding the stellar mat- 
ter onto the lens planes, although the dark-matter particles 
of the simulation represent the total matter, can then be 
compensated using extended analytic proflles with nega tive 
masses (as was done, e.g., by Wambsganss et alT]|2008 1 . 



3.4. Lensing maps and image positions 

To produce lensing maps, we set up rays starting from a 
flducial observer on a regular grid in the image plane with 
an angular fleld size and resolution suitable for the particu- 
lar application. The resulting shear and convergence maps 
may be used directly to obtain, e.g., the shear correlation 
functions and power spectra. 

We also wish to perform simulations of galaxy-galaxy 
lensing. Not only are the image positions of distant source 
galaxies affected by lensing. Also the apparent positions of 
galaxies and halos that act as lenses for background galax- 
ies (and are to be correlated with the shear fleld) are af- 
fected by lensing due to foreground matter. We therefore 
have to compute the image positions 6g given the galaxies' 

source positions (3^'' (i.e. the projected galaxy positions 
on the lens planes) and the lens mapping sampled on the 
grid of light rays in the ray-tracing algorithm. To reach 
this, we make use of a triangle technique described, e.g.. 



The ray-tracing algorithm has to compute 5 derivatives of 
the lensing potential (2 deflection angle components and 3 shear 
matrix components) starting from the matter distribution. On 
large meshes, lower-order finite differencing operations (FDs) are 
much faster than FFTs. Using FFT derivatives would require 6 
FFTs, whereas our approach requires only 2 FFTs and 5 FDs. 
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Fig. 5. The deflection angle \cy.\ as a function of the co- 
moving distance r to a point mass computed by our mesh 
method, (a) Short-range component (dashed hne) and long- 
range component (dotted line) of the deflection angle (full 
line), (b) Fractional difference between the value a calcu- 
lated by our mesh method and the theoretical value ck"^. 
For the plots, we placed 10 point masses of 8.6 x 10^h~^ Mq 
randomly on the lens plane at z — 0.5 and calculated the 
resulting deflection at 1000 random positions around each 
of them. 



Schneider et al. 



(19921. We partition the region of the 
image plane that is covered by the grid of rays into trian- 
gles formed by rays of adjacent grid points (see Fig. [t]). On 
each lens plane, we identify for each such triangle all galax- 
ies with source position inside the backtraced triangle. The 
image positions of these galaxies are then computed by lin- 
ear interpolation between the rays. This method takes into 
account strong lensing, as a galaxy might lie in more than 
one triangle on the lens plane, resulting in multiple images 
of that galaxy. 

Figure |8] illustrates how well the image positions ob- 
tained by the triangle interpolation method are mapped 
back onto the source positions by the ray-tracing. Shown is 
the difference between the true source positions and the po- 
sitions obtained by tracing back light rays starting from the 
interpolated image positions to the source plane. The differ- 
ence is always much smaller than the resolution of the mat- 
ter distribution on the lens planes. The slight anisotropy 
seen as a larger spread along one diagonal is caused by the 
particular way the image plane was partitioned into trian- 
gles. The diagonal coincides with the diagonal chosen for 
splitting the square mesh pixels into triangles. If one used 
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Fig. 6. The magniflcation for sources at z = 2 as 
a function of the angular separation to a point mass 
of lO^^/i~^M0 at redshift z = 1 computed by our mesh 
method, (a) Numerical values (symbols) compared to the 
analytical result (solid line), (b) Fractional difference be- 
tween the measured magniflcation and its theoretical 
value The magniflcation diverges at the Einstein ra- 

dius 9e — 16" (dashed vertical line). For the plots, we 
placed 10 point masses of 10^^/i~^ Mq randomly on the lens 
plane at z = 1 and calculated the resulting magniflcation 
at 1000 random positions around them. 



the other diagonal for splitting the mesh pixels into tri- 
angles, a larger spread along that diagonal would be seen 
instead. 



4. Results 

We compute various weak-lensing two-point statistics from 
a set of ray-tracing simulations and compare the results to 
semi-analytic predictions. If not stated otherwise, adaptive 
snoothing of the matter distribution on the lens planes is 
applied for the ray-tracing. 



4.1. Power spectra 

We start our discussion with the convergence power spec- 
trum Pk{^). In the flrst-order approximation ([s]) , the con- 
vergence power spectrum is given by (see, e.g., Schneider 
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image plane 



source plane 



Fig. 7. Interpolation scheme used for determining image 
positions of galaxies. The regular grid of rays in the image 
plane (left filled circles) is used to partition the image plane 
into triangles (right blue lines). The image position (left 
open circle) of a source inside a triangle (right blue lines) 
formed by the backtraced rays on the source plane (right 
filled circles) is then determined by linear interpolation. 
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Fig. 8. Comparison of the true source positions and the 
source positions obtained by tracing back light rays through 
the Millennium Simulation starting from the image po- 
sitions computed by the interpolation method. Shown is 
the difference Ax between true and traced-back comoving 
source position for 1000 galaxy centres at z = 1. At this 
redshift, 0.010ft.~^kpc comoving correspond to an angle of 
10"'^ arcsec. The right and upper frames sides are labelled 
by the corresponding angular difference between true 
and traced-back source position in units A^mosh = 1" of 
the mesh spacing used for the rays. 



20061 



dwq-'{w)Ps t{w) 



fK{w) 



Here, 



q{w) 



, , , /Jk{w -w) 
fK{w') 



(22) 



(23) 



with the probability distribution Ps{w) of visible sources 
in comoving distance. Furthermore, Ps{t,k) denotes the 
three-dimensional power spectrum of the matter contrast 
6 at cosmic time t and comoving wave vector k. For 
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Fig. 9. Convergence power spectra Pk(^) for sources at red- 
shift z — 1 (lower curves) and z — 2 (upper curves). The 
simulation results (from ~ 30 random fields of 3 x 3deg^) 
are shown as diamonds with errorbars (indicating standard 
deviation calculated from the field-to-field variance), the 
corresponding first- order predictions as solid l ines. The pre- 
dictions using the Peacock & Dodds (1996) prescription 



together with the transfer function from Eisenstein fc Hu] 
(19^ are given a s dotted lines, those obtained from the 
|Smitn et al. (20031 fitting formula as dashed lines. The pre- 
di ctions of a halomo del using the concentration parameters 
of Neto et al. ( 2007 1 are shown as dash-dotted lines. 



an accurate comparison with the results obtained by our 
ray-tracing algorithm, we use Eq. |22] together with the 
three-dimensional power spectr a Ps jk) measured fr om the 
Millennium Simulation (see also Vale fc White|2003 1. In the 
following, we will call the resulting power spectra, first- order 
prediction for brevity. 

In Fig. d we compare the first-order prediction to the 
convergence power spectra obta ined from the ra y-tra cing. 
As has already b een observed by Jain et al. ( 2000 1 and Vale 



& White ( 2003 1 , the power spectra from the ray-tracing 
are m very good agreement with the first-order prediction. 
For the considered source redshifts, the difference < 2% for 
£ < 10^. The larger deviations at wave numbers £ > 2x 10^ 
are due to smoothing effects discussed below. 

In the first-order approximation, the power spectra of 
the convergence an d the shear are ide ntical. As has already 
been found, e.g., by Jain et al. (20001, the convergence and 
shear power spectra from ray-tracmg agree very well, too. 
On scales £ > 1000, the difference between both is well 
below one percent in our ray-tracing results. 

If the first-order prediction for the convergence power- 
spectra is assumed to be correct to very high accuracy, the 
smoothing tests can be considered as a test of the accu- 
racy of our ray-tracing algorithm. Then the results shown 
in Fig. [9] suggest that the ray-tracing is able to repro- 
duce weak-lensing effects within ~ 3% accuracy on scales 
300 <^< 20 000. 

The comparison of the ray-tracing power spectra with 
some of the popular fitting f ormulae is less encourag- 



ing: Both the prescriptions by Peacock fc Dodds ( [1996 1 
(with the tra n sfer fu nction by Eisenstein & Hu 1999p and 
Smith et al.| (20031 strongly underpredict the power on 
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Fig. 10. Convergence power spectra Pk{£) for sources 
at redshift z — 1. Compared are the results from 
ray-tracing (symbols) using various smoothing schemes 
(none/Gaussian with fixed scale ^s/adaptive) and the corre- 
sponding first-order prediction (lines) obtained projecting 
and smoothing the measured 3D power spectra of the ac- 
tual mass distribution in the simulation. 



intermediate and small scales. These fitting formulae are 
based on older simulations, whose matter power spectra 
are noticeably different from the power spectra of more 
recent, higher-resolution simulations. The deviations from 
the simulated convergence power spectra exceed 30% for 
£ > 10000, so these fitting formulae seem to be of limited 
use for the interpretation of data from future weak-lensing 
surveys. 



A prediction b ased o n th e popular halo model (Seljak 



2000 Cooray & Sheth 20021) and the halo concentration- 



mass relation of Neto et aT ( 2007 1 provides a better fit to 



the convergence power spectrum. There are, however, still 
deviations (« 10%), in particular for higher source redshifts 
and intermediate scales (i.e. £ Ki 1 — 2 x 10^). This coincides 
with the transition region of the one- and two-halo terms, 
which is difficul t to model accurate ly due to halo exclusion 
et al.|2005[ and references therein) , 
in our predii 



effects (see e.g. Tinker 



which are not included 



prediction. 



As mentioned above, the deviations of the measured 
power spectra and the first-order predictions at large £ are 
due to smoothing effects. In Fig. 10 we present the conver- 
gence power spectra from ray-tracing runs of the same set 
of fields (with a cumulative area of 80 deg^ and sources at 
z = 1), but with different smoothing schemes. In addition 
to adaptive smoothing, which is intractable analytically, we 
also employ smoothing with a Gaussian kernel of fixed co- 
moving size on the lens planes. The ray-tracing simulations 
with Gaussian smoothing on the lens planes show - apart 
from sampling variance - perfect agreement with the first- 
order prediction if the smoothing is into taken into account 
there. Only the spectrum for the smallest smoothing length 
shows some aliasing effects on very small scales. The spec- 
trum of the adaptive-smoothing runs happens to match the 
spectrum for a Gaussian smoothing length of 10h~^ kpc co- 
moving quite well, but one should be cautious when consid- 
ering this as an "effective" smoothing length in a different 
context. 
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Fig. 11. Aperture mass dispersion (M|) (i?) (a) and 
(M|) {i}) (b) as a function of filter scale ■& for sources at 
z — 1 and z = 2. Compared are the first-order prediction 
(solid lines, E-mode only) and the results from ray-tracing 
(symbols with error bars indicating standard deviation, ob- 
tained from 7 fields of 5 x 5deg ). 

4.2. Aperture-mass statistics 

A suitable cosmic-shear measure that allows one to de- 
compose the shear signal in a finite-sized field into E- and 
B-modes is the aperture mass dispersion (Schneider et al. 



1998, 2002) . The E- and B-mode aperture mass at position 
& on the sky and scale are defined by: 

Ml^{9,^)^ J d^e' Q{e' -9,^) jt,x{e\e' -6). (24) 

In this w ork we use the polyno mial filter function Q pro- 



posed by Schneider et al. 



(19981: 



6|6>p 



1 - 



(25) 



The tangential and cross components of the shear are de- 
fined by 



-2i0(e) 



7x(0',0)=-3f7(e')e-''^(^' 



(26a) 
(26b) 

where 4>{6) is the polar angle for the direction defined by 
d. 
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Fig. 12. B-mode aperture mass dispersion (M|) for 
sources at z = 2 decomposed into the contributions by 
lens-lens coupling and Born corrections. Ray-tracing results 
(symbols with error bars indicating standard deviation, cal- 
culated from 7 fields of Ixldeg^): full ray-tracing (squares), 
only Born corrections (dia monds), only lens-lens coupling 
(triangles). Predictions by Cooray & Hu ( 2002[) : full sig- 
nal (dashed line) , only Born corrections (dotted line) , only 
lens-lens coupli ng (dash-dotted line). Prediction by [Hirataj 
fc Seljak|([2003|: full signal (sohd hne). 



An estimate for the aperture mass dispersion 
(M| g) {d) as a function of the filter scale d can be com- 
puted from a given shear field by a spatial average. Figure 
[TT| i shows the E-mode aperture mass dispersion measured 
from our set of simulations. The dispersion measured from 
the ray-tracing i s in very good agreem ent with the first- 
order prediction ( Schneider et al.|[l998 l: 



(M|) (t?) = 



288 



(27) 



where Pk{£) is given by Eq. (22 
function of the first kind 



and J4 denotes a Bessel 



The deviations of the measured 
aperture mass dispersion from the first-order prediction 
seen on scales < 0.5 arcmin can be attributed to smoothing. 

In the first-order approximation, the B-mode aperture 
mass dispersion {■&) vanishes. The measured B-mode 

dispersion from the full ray-tracing is shown in Fig. \TT]p . 
The B-mode signal is at least 3 orders of magnitude smaller 
than the E-mode. On larger scales their ratio even drops 
below 10~^. 

Theoretical predictions of the am plitude of the lensing- 
indu ced B-modes have been made by Cooray & Hu ( 2002 ) 
and Hirata & Seljak (2003), who calculated corrections 
to the E- and B-mode shear power spectra by expanding 
Eq. PI to second order in the gravitational potential. As 
Fig. [12] illustrates, the predictions based on their methods 
(and the measured three-dimensional power spectra of the 
Millennium Simulation) are of the correct order of mag- 
nitude and reproduce some qualitative features of the ray- 
tracing simulations, but the match is far from being perfect. 
While the B-mode predictions are lower by a factor of « 2 
on small scales, the signal measured from the ray-tracing 
declines much more quickly on larger scales. However, the 
discre pancies are not large enou gh to challenge the find- 
ing of Shapiro & Cooray (2006) that the lensing-induced 
B-mode is unimportant even for an all-sky survey. 



In order to determine their individual contributions to 
the total B-mode signal, we switch off ray deflections [i.e. 
we em ploy the Born approximation by setting 0^^) = Vfc 
in Eq. (15)] and/or lens-lens coupling [i.e. we set A*^*^^^^ = 
1 in thethird term on the r.h.s. of Eq. (16)]. Again the 
predictions and the measured signal differ-^ factors ~ 2 
on small scales, and the measured signal decreases much 
stronger with increasing scale. 

We closely examined various steps involved in the calcu- 
lations to exclude numerical artifacts as the reason for the 
discrepancy. The smoothing tests show that the smoothing 
we applied to the matter distribution on the lens planes 
can only account for deviations on scales < 0.5 arcmin. 
Examining the variance between the different ray-traced 
fields, we can exclude 'cosmic variance' as a major source 
of the discrepancy. 

The ray-tracing results did not change when different 
ways of estimating M| in real and Fourier space, as well 
as different methods of numerical integration for the theo- 
retical curves were used. Furthermore, only a tiny B-mode 
(at least 6 orders of magnitudes smaller than the E-mode) 
remained, when both ray deflections and lens-lens coupling 
were switched off in the simulation (which is essentially 
equivalent to the first-order approximation). The origin 
of this tiny signal is found to be the interpolation of the 
Jacobian matrix between the grid points to obtain their 
values at the light ray positions. Sampling a B-mode-free, 
continuous shear field on a grid and subsequent interpo- 
lation yields again a continuous shear field. This, however, 
agrees exactly with the original field only at the grid points. 
Therefore, it may in general contain a small B-mode con- 
tribution, depending on the grid resolution and the inter- 
polation scheme used. 

4.3. Galaxy-galaxy lensing 

We test the effect of Born corrections and lens-lens coupling 
on galaxy-galaxy lensing (GGL) by producing a catalogue 
of unbiased mock galaxies. We achieve this by first drawing 
a number of simulation particles on each lens plane at ran- 
dom and using their positions as lens galaxy positions in 
the algorithm described in Sec. |3.4| We then obtain a cat- 
alogue of source galaxies by randomly sampling positions 
in the image plane assuming a uniform image distribution 
over the field-of-view. 

The GGL signal we are interested in is given by the 
mean tangential shear (7t) {■&) at the image positions of the 
source galaxies as a function of angular separation to the 
positions of the lens galaxies. In the simple case of unbiased 
galaxies considered here, the expected GGL signal can be 
computed in the first-order approximation by: 



(7t) m 



1 

2tt 



Pi{w)q{w) 
fK{w) 



fK{w) 



, (28) 



where J2 is a Bessel function of the first kind, p\{'w) is the 
probability distribution of the lens galaxies' distances, the 



lensing weight q{w) is given by Eq. (23), and Ps denotes 
again the 3D matter power spectrum. For simplicity, we 
will consider a volume-limited sample of lens galaxies with 
constant comoving density in the following. 
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Fig. 13. Galaxy-galaxy-lensing signal for sources at red- 
sliift z = 1 and unbiased lens galaxies with a constant co- 
moving mean density between z = and z — I. (a) Shown 
are the measured tangential component {'-ft) (i?) of the shear 
from full ray-tracing (dia mon ds) and ray-tracing using the 
first-order appr oximation (17 1 (squares), and the first-order 
prediction ( p8| (solid line)^b) Measured cross component 
(7x)(i9) from full ray-tracing (diamonds) and first-order 
ray-tracing (squares) . Error bars denote the standard devia- 
tion calculated from a set of 24 simulated fields of 3 x 3 deg . 



Due to statistical parity invariance, the cross component 
7x is expected to vanish when averaged over many source- 
lens pairs. The observed mean cross component (7x) can 
therefore be used as a test for systematic effects and 'cosmic 
variance'. As shown in Fig. 13 (7x) is consistent with zero 
in our ray-tracing. 

While the cross component 7x provides a test for sys- 



tematic effects, the tangential shear 74 contains the desired 
information about the matter and galaxy distribution. As 
can be seen in Fig. 13 the mean tangential shear (jt) is 
significantly smaller (w 10 — 20% at an angular separation 
of 1 arcmin) in th e ra y-tracing than expected from the first- 
order prediction (28 1. 



The reason for this discrepancy is magnification bias: 
Lenses, i.e. dense matter structures such as galaxies or clus- 
ters with their dark matter halos, magnify the regions be- 
hind them. The magnification reduces the apparent num- 
ber density of higher-redshift lens galaxies around lower- 
redshift lenses in a volume limited survey (as has been 
simulated here). Underdense regions, on the other hand, 
demagnify the regions behind them, thereby increasing the 
apparent number density of lens galaxies behind them. The 



de-/magnification leads to an anticorrelation between the 
positions of high-redshift lens galaxies and the tangential 
shear induced by low-redshift structures. The anticorrela- 
tion reduces the signal (74) compared to the first-order ap- 
proximation]^ We can suppress the magnification bias in 
the ray-tracing by switching off the defiections and using 
Eq. ( 17 1 to calculate the distortions. In this case our sim- 



ulations are fully consistent with the first-order prediction, 
as is shown in Fig. [13] 

The effect of the magnification bias on the GGL de- 
pends on the redshift distribution of the sources and the 
lenses. Moreover, the shape of the lens luminosity function 
may be important if the lens population is selected using 
a magnitude limit. For example, the first-order approxima- 
tion may wnrferestimate (74) for a lens population with a 
very steep luminosity function near the survey magnitude 
limit. We reserve a more detailed investigation of this effect 
with realistic source and lens distributions for future work. 



5. Summary 

In this work, we have described a new variant of the 
multiple-lens-plane algorithm, which is particularly suited 
for ray-tracing through very large cosmological A^-body 
simulations. The algorithm differs in some important de- 
tails from previous works. This allows us to take full advan- 
tage of the unprecedented statistical power offered by the 
large volume and high spatial and mass resolution of the 
Millennium Simulation. The features discussed include: a 
tilted line-of-sight (to avoid periodic repetition of structures 
along the line-of-sight), adaptive slice boundaries (to avoid 
the slicing and duplication of bound structures), adaptive 
smoothing of the projected matter distribution on the lens 
planes (to reduce shot noise from the particles) , a mutliple- 
mesh method for calculating the light deflections and dis- 
tortions at the lens planes (which takes into account the 
small-scale and large-scale structure simultaneously), and 
a method to include galaxies (as lenses and sources) from 
semi-analytic galaxy-formation models in the ray-tracing 
process. 

We have used the ray-tracing code and the Millennium 
Simulation to investigate the impact of lens-lens coupling 
and multiple ray deflections on various cosmic shear two- 
point statistics. We have computed convergence power 
spectra from a set of ray-tracing realisations. For testing 
and comparison, we have also computed a first-order predic- 
tion of the convergence power spectrum using the measured 
three-dimensional power spectra of the mass distribution 
in the Millennium Simulation. We find that this first-order 
prediction agrees very well with the ray-tracing results ex- 
cept for very small scales (the difference is > 5% only for 
£ > 20000), where smoothing on the lens planes becomes 
important. 

Comparing the convergence power spectrum from the 
ray-tracing to the predictions based on the fitting formu- 
l ae for the matter power spec trum by [Peacock fc Dodds| 
(1996) and Smith et al. ( 2003[ ), we find significant discrep- 
ancies (> 30% for I > 10000), casting the usefulness of 
these fitting formulae for cosmological parameter estima- 



Note that in the first-order approximation, magnification ef- 
fects are neglected. Thus, the positions of galaxies at any given 
redshift are uncorrelated with the shear induced by galaxies at 
different redshifts. 
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tion for future surveys into doubt. A prediction based on 
the popular halo model and the halo concentration-mass 
relation of Neto et al. (2007) fits better, but there are still 
noticeable deviations, in particular for higher source red- 
shifts 10% for sources at = 2). This indicates a need 
for more accurate descriptions of matter power spectra. 

Furthermore, we have computed the E- and B-mode 
aperture mass dispersion using our ray-tracing algorithm. 
We find the B-mode to be finite, but at least three orders 
of magnitude smaller than the E-mode. The amplitude of 
the B-mode is slightly larger and shows a different scal e de- 
pendence than the pre dictions of Cooray & Hu ( 2002 1 and 
Hirata fc Seljak ( 2003 1 . We have performed various tests to 
exclude numerical artifacts as the origin of the deviations. 
Despite these discrepanc ies, we can confirm the finding of 
Shapiro fc Cooray ( 2006 1 that the lensing-induced B-mode 
can be safely neglected even in an all-sky survey. 

Corrections to the first-order approximation can have a 
considerable impact on galaxy-galaxy lensing. In the simple 
case of a volume-limited sample of unbiased lens galaxies 
and all sources at redshifts z = 1, the first-order approx- 
imation overestimates the mean tangential shear around 
lenses by ~ 10 — 20% at an angular separation of 1 arcmin 
due to its failure to incorporate the magnification bias. The 
impact of the magnification bias on the galaxy-galaxy lens- 
ing signal depends on the survey selection criteria and the 
luminosity and redshift distribution of the sources and the 
lenses. A detailed investigation of this effect should be car- 
ried out in future work. 
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Appendix A: Lattice planes 

The periodicity of our matter distribution along and per- 
pendicular to the line-of-sight ca n be studied within the the - 
ory of crystal lattices (see, e.g., Ashcroft fc Mermin||1976 1. 
Here, we give a practical explanation rather than a rigorous 
proof. 

Consider an array of unit cubes forming a simple cu- 
bic lattice with lattice constant unity. Choose two linearly 
independent lattice vectors p and q with p — {pi,P2,Pi) 
and q = ((7i, (72, 93) and Pi,qi £ 1>. These two vectors span 
a plane which is perpendicular to the lattice vector n with 
n = (711,712,713) = {pi,P2,P3) X (91,92,93), n-i e Z. 

Since the plane-spanning vectors are lattice vectors, the 
plane is itself periodic and therefore represents a plane lat- 
tice. With p and q as basis vectors of the plane lattice, the 
plane is periodic along the direction of p and q with peri- 
odicity length IpI and \q\, respectively. The parallelogram 
constructed from p and q represents a unit cell of the plane 
lattice with a cell area of |p x q| = \n\. One can show that 
there is no smaller unit cell if the integer coefficients rti, 
712, and 713 are coprime. Furthermore, there is no shorter 
non-zero lattice vector perpendicular to the plane than n 
in this case, and hence, the shortest periodicity along the 
normal direction is \n\. 

For the computational cube of the Millennium 
Simulation with side length L — 500h~^ Mpc, the lengths 



and areas above have to be multiplied by L and L^, re- 
spectively. Our choices p = L(3, —1, 0) and q = L{1, 3, — 1) 
yield a LOS vector n = L{1, 3, 10) with |77| = 5.244/i-i Gpc 
and a rectangular area of 1.581/i~^ Gpc x 1.658/i~^ Gpc for 
the lens planes. 
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